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Abstract 

In order to address the imprecision often introduced by widening operators in static analy¬ 
sis, policy iteration based on min-computations amounts to considering the characterization of 
reachable value set of a program as an iterative computation of policies, starting from a post- 
fixpoint. Computing each policy and the associated invariant relies on a sequence of numerical 
optimizations. While the early research efforts relied on linear programming (LP) to address 
linear properties of linear programs, the current state of the art is still limited to the analysis of 
linear programs with at most quadratic invariants, relying on semidefinite programming (SDP) 
solvers to compute policies, and LP solvers to refine invariants. 

We propose here to extend the class of programs considered through the use of Sums-of- 
Squares (SOS) based optimization. Our approach enables the precise analysis of switched sys¬ 
tems with polynomial updates and guards. The analysis presented has been implemented in 
Matlab and applied on existing programs coming from the system control literature, improving 
both the range of analyzable systems and the precision of previously handled ones. 


1 Introduction 

A wide set of critical systems software including controller systems, rely on numerical computations. 
Those systems range from aircraft controllers, car engine control, anti-collision systems for aircrafts 
or UAVs, to nuclear powerplant monitors and medical devices such a pacemakers or insulin pumps. 

In all cases, the software part implements the execution of an endless loop that reads the sensor 
inputs, updates its internal states and controls actuators. However the analysis of such software is 
hardly feasible with classical static analysis tools based on linear abstractions. In fact, according 
to early results in control theory from Lyapunov in the 19th century, a linear system is defined 

*The work was done when the author was at IRIT, Universite Paul Sabatier at Toulouse and was supported by 
the CIMI (Centre International de Mathematiques et d’Informatique) Excellence program ANR-11-LABX-0040-CIMI 
within the program ANR-ll-IDEX-0002-02 during a postdoctoral fellowship. 
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as asymptotically stable iff it satisfies the Lyapunov criterion, i.e. the existence of a quadratic 
invariant. In this view, it is important to develop new analysis tools able to support quadratic or 
richer polynomial invariants. 
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(a) A one-loop program with a switch-case loop body. 


Let X° = {y | ri(y) M 0 and ... and r^ Q {y) M 0} 
and for all cases i, X 1 = {y \ r\{y) M 0 and ... and r l n (y) M 0}. 

We can define the discrete-time switched system: 

x 0 G X in , VfceN, x k+1 = T(x fc ), where T(x) = T\x) if x G X i n X° 
(b) The discrete-time switched system correspondence of the program. 


Figure 1: One-loop programs with switch-case loop body and its representation as a switched 
system. 

While most controllers are linear, their actual implementation is always more complex, e.g. in 
order to cope with safety, additional constructs such as antiwindups or saturations are introduced. 
Another classical approach is the use of linear parameter varying controllers (LPV) were the gains 
of the linear controller vary depending on conditions: this becomes piecewise polynomial controllers 
at the implementation level. 

We are interested here in bounding the set of reachable values of such controllers using sound 
analyses, that is computing a sound over approximation of reachable states. We specifically focus 
on a class of programs larger than linear systems: constrained piecewise polynomial systems. 

A loop is performed while a conjunction of polynomial inequalities is satisfied. Within this loop, 
different polynomial updates are performed depending on conjunctions of polynomial constraints. 
This class of programs is represented in Fig. la The program can be analyzed through its switched 
system representation (see Fig 


lb) 


In the obtained system, the conditions to switch are only 
governed by the state variable: at each time k, we consider the dynamics T l such that x k G X 1 . 
So, the switch conditions do not depend on the time (we do not consider the time when we reach 
X*) and are not defined from random variables. 

From the point-of-view of the dynamical systems or control theory, to find an over-approximation 


of reachable values set of a program of the class described in Fig. la can be reduced to find a positive 


invariant of its switched system representation (see Fig. lb). 


Moreover, the class of switched systems where the switch conditions only depend on the state 
variable is classical in control theory and includes the nonlinear systems with dead-zone, saturation, 


is either the strict < or the weak (<) comparison operator over reals. 
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resets or hysteresis. 


Related Works 

Reachability analysis is a long-standing problem in dynamical systems theory, especially when the 
system dynamics is nonlinear. In the particular case of polynomial systems, this problem has 
recently attracted several research efforts. 

In |WLW13j . the authors use a method based on sublevel sets of polynomials to analyze the 
reachable set of a continuous time polynomial system with initial/general constraints being en¬ 
coded by semialgebraic sets. Their method relies on a so-called iterative advection algorithm 
based on Sums-of-Squares (SOS) and semidefinite programming (SDP) to compute either inner 
of outer approximation of the backward reachable set, also named domain of attraction (DoA). 
The stability analysis of continuous-time hybrid systems with SOS certificates was investigated 
in | jPP09) . (PTT1 fi| have recently applied analogous techniques to perform stability analysis and 
controller synthesis in the context of robotics. Other studies rely on SOS reinforcement and moment 
relaxations to obtain hierarchies of approximations converging to sets of interest such as the DoA 
in continuous-time, either from outside in | |HK14| or from inside in [KHJ12] , This approach relies 
on a linearization of the ordinary differential equation involved in the polynomial system, based 
on Liouville’s Equation satisfied by adequate occupation measures. This framework has been ex¬ 
tended to hybrid systems in [SVBT14] , as well as to synthesis of feedback controllers in }MVTT13 . 
Liouville’s Equation can also be used to approximate other sets of interest, such as the maximal 
controlled invariant for either discrete- or continuous-time systems (see }KHJ14j ). 

By contrast with these prior works, our approach focuses on computing over approximation of 
the forward reachable set of a discrete-time polynomial system with finitely many guards, thanks 
to an algorithm relying on SOS and template policy iterations. 

Template abstractions were introduced in }SSM04| as a way to define an abstraction based 
on an a-priori known vector of templates, i.e. numerical expressions over the program variables. 
An abstract element is then defined as a vector of reals defining bounds bi over the templates pp. 
Pi(x i,...,x d ) < bi . 

Initially templates were used in the classical abstract interpretation setting to compute Kleene 
fixpoints with linear functions pi- Typically, the values of the bound bi increases during the fixpoint 
computation until convergence to a post-fixpoint. 

Later in AGO. 10 the authors proposed to consider generalized templates including quadratic 
forms and compute directly the fixpoint of these template-based abstractions using numerical opti¬ 
mization. When considering the sub-class of linear programs, fixpoints can be computed by using 
a finite sequence of linear (LP) and semidefinite (SDP) optimization problems. 

Two dual approaches, respectively called Max-policies and Min-policies, can be applied. Max- 
policies |GS07| iterate from initial states and compute policies as relaxations through rewriting of 
an optimization problem (forgetting about rank conditions). Min-policies jGGTZ07l lAGGlO] rely 
on duality principle and determine a policy through the computation of a Lagrange multiplier. 

Contributions 

The present paper is a followup of |AGM15j . in which Sums-of-Squares (SOS) programming is 
used to analyze properties (such as boundedness or safety) of piecewise discrete polynomial sys¬ 
tems. The main contribution is the extension of the Min-policy iteration algorithm to improve the 
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precision of the analysis of boundedness for such systems. Here, we develop a policy iteration on 
template domains based on polynomials. The current approach improves the previous frameworks 
developed in [GGTZ071 IAGG12| to handle affine (or very specific) piecewise affine systems with 
quadratic templates and semidefinite programming. This improvement comes from the use of SOS 
programming to develop an SOS extension of a relaxed functional, sharing the same properties as 
the one defined in [AGG12j . In particular, we show, as in IMT41 , that this policy iteration has 
the desired convergence guarantees. 

Organization of the paper 

The paper is structured as follows: we characterize the class of programs considered - constrained 
piecewise discrete-time polynomial systems - together with their collecting semantics as a least 
fixpoint. Then, in Section [3j we recall definitions of generalized templates, their expression as 
an abstract domain and the definition of the abstract transfer function. Section [4] proposes an 
abstraction of the transfer function using an SOS reinforcement, while Section [5] relies on this 
abstraction to perform policy iteration. Finally Section [6] presents experiments. 

2 Constrained Piecewise Discrete-Time Polynomial Systems: Def¬ 
inition and Collecting Semantics 

In this paper, we are interested in proving automatically that the set of values taken by the vari¬ 
ables of the analyzed program is bounded. In the rest of the paper, we analyze the program by 
decomposing it using its dynamic system representation. The boundedness problem is thus reduced 
to prove that the set of all the possible trajectories of a dynamical system is bounded. Since the 
analyzed program has the form depicted by Figure [laj we consider the special class of discrete-time 
dynamical systems introduced in Figure [Tb| that is: 

(i) their dynamic law T is a piecewise polynomial function, and 

(ii) the state variable x is constrained to live in some given basic semi-algebraic set0 

We recall that a set is a basic semi-algebraic set if and only if it can be represented as a 
conjunction of strict or weak polynomial inequalities (“basic” means that no disjunction occurs). 
Note that T is a piecewise polynomial function with respect to a given partition, meaning that if 
we restrict T to be an element of the partition then T is a polynomial function. 

We now give a formal definition of constrained piecewise discrete-time polynomial system (PPS 
for short). 

First to define a PPS, we need a partition. Let X be a finite set of partition labels and X = 
{X 1 C | i £ X} be a partitioning, that is a given family of basic semi-algebraic sets satisfying 
the following: 

[Jr = R d , Vi, j G X, i ± j => X i n X j = 0 . ( 1 ) 

i£l 

Each set X 1 of the partition corresponds to a case in the loop body of the program given in Figure [la| 
as the cases are assumed to be disjoint. By definition of basic semi-algebraic sets, it follows that 

J For instance membership of the sub-level set {x £ R d | 1 — 11*111 < 0}, thus this does not entail boundedness of 
variable values. 
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for all i 6l, there exists a family of n* polynomials {r* ; j E [n*]} such that: 

X i = {x E R d | rj(x) MOVjE [n*]} . (2) 

where N is either < or < and [n*] denotes the set {1 , Eq. © matches with the program 
given in Figure la Guards are defined as conjunctions of polynomial inequalities and thus are 
basic semi-algebraic sets. 

The second tool needed to define a PPS is the piecewise polynomial dynamic relative to the 
partition. Let T : i —> be a piecewise polynomial function w.r.t. to the partitioning X. By 

definition, there exists a family of polynomials {T* : H > M rf , i E 1} such that for all iEl: 


T(x) = T(x), Vx E X 
Eq. ([3]) matches with the polynomial updates in Figure la 


(3) 


Finally, it remains to define the initialization and the set where the state variable lies. Let 
X m and A" 0 be two basic senri-algebraic sets of A m supposed to be compact, i.e. closed and 
bounded. The two sets can be represented as in Eq. ([ 2 ]) using their respective family of ni n and no 
polynomials: 


X m = iE 


jx E R d | r) n (x) N 0 V j E [ni n ]| 


and A 0 = |iE 


jx E R d | r°(x) N 0 V j E [no] j , 


where for all j E [n; n ] , rfi : 




and for all k E [no], r® : M“ 1 —> M is a polynomial. 


The set X m and X u respectively denote the set of initial states of the program and the set 


which defines the loop condition in Figure la 


Let X be the family of sets {X\i E X} satisfying Eq. 0 and T be the family of functions 
{T\i E X} satisfying Eq. Q. We define the PPS associated to the quadruple (X m , A 0 , X,T) as 
the system satisfying the following dynamic: 

xo E X in , and V k E N, if Xfc E A 0 , Xk+i = X(x&). (4) 

In the rest of the paper, a PPS dynamical system is identified by a quadruple (X m , A 0 , X , T). 

Example 2.1 (Running example). We consider the following running example corresponding to a 
slightly modified version of IAJ131 Example 3]. By comparison with 1 AJ13 . Example 3], the semi- 
algebraic set A 1 (resp. X 2 ) is introduced to represent conditions under which we use the polynomial 
update T 1 (resp. T 2 ). The PPS is the quadruple (X in , X°, {X 1 , X 2 },{T 1 ,T 2 }), where: 


X in = [—1,1] x [—1,1] 

A 0 = M 2 


and 


A 1 = {x E 
A 2 ={xE 


| -x\ + 1 < 0} 
£ 2 j x\ - 1 < 0} 


and the family of functions { T 1 ,T 2 }, defined as follows: 


. / / 0.687xi + 0.558x2 — 0.0001x1X2^ 

X i (x 1 , x 2 ) = | nono _ , n I 


-0.292x1 + 0.773x2 
and 


_ 2 / \ ( 0.369xi + 0.532x9 — O.OOOlx? \ 

T (X1 ’ X2)= ^-1.27xi+ 0.12X2-O.OOOIX 1 X 2 J 


Its (discretized) reachable value set is simulated and depicted at Figure [I| 
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Figure 2: Running example simulation 


We recall that our objective is to prove automatically that the set of the 
of the system is bounded. This set, also called the reachable values set or the 
of the system, is defined as follows: 

9T = (J T^ o (X in ) 

fee N A 


possible trajectories 
collecting semantics 

(5) 


n (k) 

where XI ,, is the restriction of T over X u .The notation X v 

1x0 | x 0 


stands for composing /c-times the map 


T \x°' 

To prove this boundedness property, we can compute and do some analysis to prove that iR 
is bounded. Nevertheless, the computation of iR cannot be done in general and instead, we have to 
compute an over-approximation of iR and show that this approximation is bounded. 

The usual abstract interpretation methodology to characterize and to construct this over¬ 
approximation relies on the representation of 1R as the smallest fixed point of a monotone map 
over a complete lattice. In other words, iR satisfies: 


m = T(9d nx°) u x' m = (J T*(iR ni°n x*) ux in . 

i£l 


Let us define p(M d ) the set of subsets of and introduce the function F : p(M rf ) i—?> p(M rf ) defined 
for all C £ p(M d ) by: 

F(C) = T(c nx°) ux in = (J t\c m°n x { ) ux in . (6) 

i£l 

Thus, iR is the smallest fixed point of F and from Tarski’s Theorem, since F is monotone and 
p(M rf ) is a complete lattice: 


m = min{C £ p(M d ) | F(C) C C} . 


(7) 
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Finally to compute an over-approximation of 9a it suffices to compute a set C such that F{C ) C C. 
A set C which satisfies F(C) C C is called an inductive invariant 

The rest of the paper addresses the computation of a sound over-approximation of 91 using its 
definition as the smallest fixpoint of F (Eq. ©)■ 


3 Basic Semi-algebraic Inductive Invariants Set 


An easy way to over-approximate the set of reachable values is to restrict the set of inductive 
invariants that we consider. We propose to restrict the class of such invariants to basic serni- 
algebraic sets using template abstractions. A template abstraction consists of representing a given 
set as the intersection of sublevel sets of a-priori fixed functions depending on the state variables. 
Such functions are called templates. Then computing an inductive invariant in the templates 
domain boils down to providing, for each template p, a bound w(p) such that the intersection over 
the templates p of sublevel sets {x G | p(x) < w(p)} is an inductive invariant. In our context, a 
template is simply an a-priori fixed multivariate polynomial. 

The overall method is not new and corresponds to a specialization of the templates abstraction 
(see jAGGlOl 1AGG12] ) to polynomial templates. However, in practice, the method developed 
in |AGG1Q[ IAGG12j is restricted to template polynomials of degree 2 (quadratic forms) and affine 
systems or a very restricted class of piecewise affine systems. 

Next we give formal details about the polynomial template abstraction and the equations that 
must satisfy the template bounds vector w to generate an inductive invariant. From now on, we 
denote by P the set of templates and by F (V, the set of functions from P to M = MU{— oo, +oo}. 

We equip F (V,M^ with the functional partial order <jr i.e. v <T w iff v(p) < w(p) for all p £ P. 
Let w G F (V, . The sets that we consider are of the form: 

w* = {x € | p(x) < w(p),Vp G P} . (8) 


Example 3.1. Let us define qi(x) = qi(xi,X 2 ) = xf, q 2 (x) = q 2 (xi, X 2 ) = x\ and let us con¬ 
sider a well-chosen polynomial p of degree 6. We will explain in Subsection \5.3\ how to generate 
automatically this template p. Let us define P mn := {qi,q 2 ,p}- 


Consider the function w u over 


w°(qi) = w°(q 2 ) = 2.1391 and w a (p) = 0, the set w° = 


{(xi, x 2 ) G M 2 | xf < w°(qi), x '2 < w°(q 2 ), p(x) < w°(p)} is presented in Figure [3j 

Now let us take the function w 1 over P run defined by w 1 (qi) = 1.5503, w l {q 2 ) = 1.9501 and 
w 1 (p) = 0, the set w 1 * = {(xi,X 2 ) G M 2 | x\ < w 1 (qi), x\ < ud(g 2 ), p(x) < w 1 (p)} is also 
presented in Figure [5[ 

Since the set of black dots in Figure [?] belongs to w°* and w 1 *, we guess that w°* and w 1 * 
contain both the reachability set of the system from Example\2. 1\ To formally prove it, one way is 


to show that they are also inductive invariants for this system. 


We restrict the class of inductive invariants to those of the form (|8]) and characterize the 
inductiveness for such sets. Since each polynomial template p is fixed, the considered variables 
that we handle are the template bounds w G F (V, . Therefore, we need to translate the 

3 In the dynamical systems theory, the inductive invariant sets are called positive invariant. 
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The semialgebraic set of gray denotes the set of Example 3.1 since p(x) < 0 is included in 


xf < 2.1391, for i E [1,2]. The dark gray region denotes the semialgebraic set w° . Black dots are 
actual reachable states of obtained by simulation. 

Figure 3: Semialgebraic sets w* for Example 


3.1 


inductiveness of the sets w* into inequalities on w. By definition the set w* is an inductive invariant 
iff F(w*) C w*, that is: 

U T\w* ni°n x*) ux in c w* . 

i£l 

By definition, w* is an inductive invariant iff: 

y p e p, Vx e (J T(w* m°n x*) ux in , P (x) < w( P ) . 

iel 

Using the definition of the supremum, w* is an inductive invariant iff: 

VpE P, sup p(x)<w(p) . 

ie[jTVn x° n x l ) u x in 

i£l 

Now, let us consider p E P. Using the fact that for all A,B C and for all functions /, 

sup / = sup{sup /, sup /}: 

AUB A B 

sup p(x) = sup < sup sup p{x), sup p(x) 

x^{j T i {w*nx°nx i )vx' m {x&ryw *nxonxq xex in 

By definition of the image: 

sup p(x) = sup < sup sup p{T l {y)), sup p(x) 

x& U,exT i (w*nx°nx i )ux in y£w* nx°nx i iex in 
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Next, we introduce the following notation, for all F: 

Ff(w)(p) := sup p(T l (x )) and X m \p) := sup p(x) . 

xe«)*nx i n.Y 0 xex in 

Finally, we define the function from F ^P, to itself, for all w E F ^P,R): 


F\w) := sup < sup pf (w), X 1: 


iex 


Note that ^ 1 correspond exactly to the notations used in |AGG10| . By construction we obtain the 
following proposition: 


Proposition 3.1. Let w £ F (P, MJ. Then w* is an inductive invariant (i.e. F(w*) C w*) iff 
F$(w) w. 


From Prop. 


3.1 


w* of the form ([8]). 


infjrc € F (P,M^ | F^{w) <jr re} identifies the smallest inductive invariant 


3.1 


i.e. 


Example 3.2 

templates basis P run from Example 
is a well-chosen polynomial of degree 6. Let w £ 
have: 


Let us consider the system defined at Example \2.1\ Let us consider the same 

= {qi,q 2 ,p} where qi(x) = x'\, q-fix) = x\ and p 
F ^P run ,Mj. For i = 1 and the templates q\, we 


F\ {w){qi) = sup (0.687xi + 0.558x2 — O.OOOIX 1 X 2) 2 . 

-xl+l<0 

X\<w{qL), x\<w{q2), p(x)<w{p) 

Indeed, X 1 = {x E M 2 | — x\ + 1 < 0} and X° = M 2 and the dynamics associated with X 1 
is the polynomial function T 1 defined for all x E M 2 by: T 1 (x) = ( 0 ' 687x fio292x^+0 7^3x2 XlX2 \ 
Since q\ computes the square of the first coordinates, this yields q\(T 1 {x)) = (0.687xi + 0.558x2 — 
O.OOOIX 1 X 2 ) 2 . 


With w E P 7 (P, Mj, computing F\w) boils down to solving a finite number of nonconvex 
polynomial optimization problems. General methods do not exist to solve such problems. In 
Section [ 4 J we propose a method based on Sums-of-Squares (SOS) to over-approximate F^{w). 


4 SOS-based Relaxed Semantics 

In this section, we introduce the relaxed functional on which we will compute a fixpoint, yielding 
a further over-approximation of the set 91 of reachable values. This relaxed functional is con¬ 
structed from a Lagrange relaxation of maximization problems involved in the evaluation of F$ and 
Sums-of-Squares strengthening of polynomial nonnegativity constraints. First, we recall manda¬ 
tory background related to Sums-of-Squares and their application in polynomial optimization. The 
interested reader is referred to }Lasf)9j for more details. 
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4.1 Sums-of-Squares Programming 

Let R[x] 2m stands for the set of polynomials of degree at most 2m, and S[x] C M[x] be the cone of 
Sums-of-Squares (SOS) polynomials, that is S[x] := { Yhi ( i%i with qi E M[x] }. 

Our work will use the simple fact that for all q E £[ 2 ], then q{x) > 0 for all x E R d as the set 
£[ 2 ] contains only nonnegative polynomials. In other words, for any given polynomial q , we can 
strengthen the constraint of q being nonnegative into the existence of an SOS decomposition of q. 

For q E M[x] 2 m , finding an SOS decomposition = Q valid over R d is equivalent to solve 

the following matrix linear feasibility problem: 

q(x) = bm(x) T Q b m (x) , Vx E R d , (9) 

where b m (x) := (1, x\,... , Xd, 2 2 , x\x 2 , ■ ■ ■, x™) (the vector of all monomials in x up to degree m) 
and the Gram matrix Q, being a semidefinite positive matrix (i.e. all the eigenvalues of Q are 
nonnegative). The size of Q (as well as the length of b m ) is ( d+ J n )- 

Example 4.1. Consider the bi-variate polynomial q(x) := 1 + x\ — 2x\x 2 + x\. With b\{x) = 
( 1 , 21 , 22 ); one looks for a semidefinite positive matrix Q such that the polynomial equality q( 2 ) = 
b\(x) T Q 61 ( 2 ) holds for all 2 E M 2 . The matrix 

10 0 \ 

0 1 -1 
0-1 1 J 



satisfies this equality and has three nonnegative eigenvalues, which are 0, 1, and 2, respectively asso¬ 
ciated to the three eigenvectors eo := (0, l/\/2, l/\/2) t, e\ := (1,0,0)t ande 2 := (0, l/x/2, —l/\/2) T . 


Defining the matrices L := (ei e 2 eo) = 




and D = 


f 0 2 0 V one obtains the decompo- 
Vo 0 0 / 


sition Q = L J D L and the equality q( 2 ) = (L 61 ( 2 )) T D (L 61 ( 2 )) = 17 ( 2 ) = 1 + (21 — 22 ) 2 , for all 
x E M 2 . The polynomial a is called an SOS certificate and guarantees that q is nonnegative. 


In practice, one can solve the general problem ([9]) by using semidefinite programming (SDP) 
solvers (e.g. Mosek [ AAOO I. SDPA |YFN + 10j . CSDP |Bor99] ). For more details about SDP, we 
refer the interested reader to |VB94j . 

The SOS reinforcement of polynomial optimization problems consists of restricting polynomial 
nonnegativity to being an element of £[ 2 ]. In case of polynomial maximization problems, the SOS 
reinforcement boils down to computing an upper bound of the real optimal value. For example 
let p E M[x] and consider the unconstrained polynomial maximization problem sup{p( 2 ), 2 E M rf }. 
Applying SOS reinforcement, we obtain: 


sup{p( 2 ), 2 E R d } = infjry | rj — p(x) > 0} < inf{?y | p — p{ 2 ) E £[ 2 ]} . (10) 


Now, let p,q E M[ 2 ] and consider the constrained polynomial maximization problem: sup{p( 2 ) \q(x) < 
0, 2 E M d }. Let A E 5H[ 2 ], then: 


sup p( 2 ) < sup p(x) — X(x) • q( 2 ) . 

q(x)<0, xSM 12 xeK d 

Indeed, suppose q( 2 ) < 0, then —X(x)q(x) > 0 and p( 2 ) < p( 2 ) — X(x)q(x). Finally taking the 
supremum over {2 E R d \ q( 2 ) < 0} provides the above inequality. Since sup{p( 2 ) — A( 2 ) • q{x), x E 
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M rf } is an unconstrained polynomial maximization problem then we apply an SOS reinforcement 
(as in Eq. ©) and we obtain: 

sup p(x) < sup p(x) — X(x) • q(x) < inf {77 | p — p — Xq E E[x]} . 

q(x)<0, a;SR d a;SM d 


Finally, note that this latter inequality is valid whatever A 6 E[i] and so we can take the infimum 
over A £ £[&] which leads to: 


sup 

q(x)<0, a:eR d 


p(x) < inf 
AeE[x] 


sup p(x) — X(x) 


q(x) < inf 

T)—p— AgeE[x] 
AeE[x] 


V- 


( 11 ) 


In Eq. (11), A is an SOS polynomial but to exploit linear programming solvers in policy iterations 


(see the fourth assertion of Prop. |5.2[ ) we restrict A to be a nonnegative scalar and in this case, since 
positive scalars are sum-of-squares polynomials of degree 0 , we obtain a safe over-approximation of 
the right-hand-side of Eq. ©• 

In presence of several constraints, we assign to each constraint an element a £ £[x], and we 
consider the product of a with its associated constraint and then the sum of all such products. 
This sum is finally added to the objective function. 

The use of such SOS polynomials for constrained polynomial optimization problem can be seen 
as a generalization of the S-Drocedure from |Yak77|. We refer to [Las m or |Par03l for applications 
in control. Note that the existence of SOS decompositions of positive polynomials over compact 
sets is ensured by the Putinar Positivstellensatz from |Put93| . 


4.2 Relaxed semantics 

The computation of F$ as a polynomial maximization problem cannot be directly performed us¬ 
ing numerical solvers. We use the SOS reinforcement mechanisms described above to relax the 
computation and characterize an abstraction of FK 

We still assume the knowledge of the template basis P, involving polynomials of degree at most 
2m. Let us define T (P,M+) the set of nonnegative functions over P i.e. g £ T (P, M_|_) iff for all 
p £ P, g{p) £ M+. Let p £ P and w £ F (P, Mj. Starting from the definition of F'. one obtains the 
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following: 

( F i( w ))(p) = 


< 


< 


sup p{T l {x)) 

q(x)<w(q), VgSP 
r}(x)<0, VjeK] 
r ° 0)5:0, Vfce[no] 

, inf , SU P P{T\x)) + J2 x (q){™(q) - ?(®)) 

Ae>(IP,K+) a;SR d . n(= ro n 

tre5[x],weE[a;],7 ;eS[x] n 'l n 

deg(oj<2mdeg T* 
deg (/zj r j )<2m deg T* 
deg( 7 ir ; °)< 2 m degT* 

inf r/ 

i,v 


-Y m( x Yi( x ) - 55 7/0*0 rV) 


1=1 


1=1 


V - p o T l - Y Kq)(w(q) -q) + Y w r l + 51 7 l r i = V , 

qSP 1=1 1=1 

S-t ' * AS J (P,M + ), cr E E[x], pi E £[x], 7/ E £[x], i)El, 
deg(<r) < 2m deg T*, 

. deg (pir\) < 2m deg T*, deg^rf) < 2mdeg T*. 

(using an SOS reinforcement to remove the sup) 

We denote by S[x] n the set of n-tuples of SOS polynomials. For clarity purpose, the dependency 
on i is omitted within the notations of the multipliers pi and 7 p Moreover, let us write Ya= 1 
(resp. )CJLi H r T) as {pi r% ) (resp. ( 7 , r 0 }). Finally, we write ( Fj l {w)){p ) the over-approximation of 
(Ff(w))(p), defined as follows: 

(F^-(w)){p) = inf p 

A,cr,^,7,77 

V - p O T l - Y Kq)( w (q) -q) + (v, r l ) + ( 7 , r°) = a 

qSP 

AsJ(P,i + ), <tES[x], p£Z[x] ni , 7ES[i]"«, p E R, 
deg(cr) < 2mdegT®, deg((/x, r l ) + ( 7 , r 0 )) < 2m deg T*. 


s. t. 


( 12 ) 


In Equation (12), the notation A is a vector of Lagrange multipliers. Each multiplier is associated 


with a constraint constructed from a template i.e. a constraint q{x) — w(q) < 0. We also introduce 
the vector of SOS polynomials p and 7 . Their role is to take into account the presence of the 
constraints x E X i and x E X° in the computation of (F^'(w))(p). Recall that X 1 and A 0 are 
basic semi-algebraic sets, then the size of the vectors p and 7 are equal to the number of polynomials 
defining X 1 and X°. 

We conclude that, for all i E I, the evaluation of Fj^ can be done using SOS programming, 
since it is reduced to solve a minimization problem with a linear objective function and linear 
combination of polynomials constrained to be sum-of-squares. 

Note that Fj^ defined at Eq. (12) is the SOS extension of the relaxed function defined in (AGG12] . 
Indeed, considering the special case where T l is affine, the templates p. q and the test functions r l , 
r° are quadratic, the vectors p i and 7 * are restricted to be nonnegative scalars, then Fj^ corresponds 
to the relaxed function defined in |AGG12l at Eq. (3.12). 


Example 4.2. We still consider the running example defined at Example \2.1\ and take again the 
same templates basis P, 


of Example 3.1 composed of q\ : 1 E ij and q 2 


a: E rj' and a 
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well-chosen polynomial p of degree 6. For the index of the partition i = 1. Recall that T l (x) = 
(°.6 8 7x^+0^558x2—o^ooo 1 x\x. 2 ) and x 1 = {x G M 2 | -s? + 1 < 0} and thus r\{x) = -xj + 1. Let 

w G F ^P run ,M^, then: 

{F^{w)){q i) = 

inf rj 

A, <7, fl,T] 

{ rj — (0.687xi + 0.558x2 — O.OOOIX 1 X 2) 2 — X(qi)(w(qi) — x 2 ) 

-\(q 2 )(w(q 2 ) -x§) - X(p)(w{p) -p(x)) +/i(x)(l - xf) = <r(x) 

AG J(P,R+), <r G E[x], fi £ S[x], ryGM, 
deg(cr) < 6 , deg(//) < 6 . 

/n practice, one cannot find any feasible solution of degree less than 6, thus we replace the degree 
constraint by the more restrictive one: deg(cr) < 6 , deg(/ 1 ) < 6 . 

The computation of F N requires the approximation of X 111 ^ := sup{p(x),x G X in }. Since X m is 
a basic semi-algebraic set and each template p is a polynomial, then the evaluation of A 111 ^ boils 
down to solving a polynomial maximization problem. Next, we use SOS reinforcement described 

. f .7^ 

above to over-approximate X m with the set X m , defined as follows: 


X' mn (p) := inf < rj 


7 1 ~ P + (v nin ,r niu ) = (To, 
i)£l, uo G E[x], i/ in G E[x] Tlin , 
deg(cro) < 2 m,deg((i/ nin , r niu )) < 2m 


n 7^. / 


Thus, the value of A in,v (p) is obtained by solving an SOS optimization problem. Since X m is a 
nonempty compact basic semi-algebraic set, this problem has a feasible solution (see the proof of 


|Las011 Th. 4.2]), ensuring that A 111 ^^) i s finite valued. 


Example 4.3. The initialization set X m of Example 2.1 is [—1,1] x [—1, 1]. It can be written as: 
{(xi, X 2 ) G M 2 | x 2 — 1 < 0, x 2 — 1 < 0}. Then, considering the same template basis of Example f.2 
and the template q\: 


IK 


A m ' v ( gi ) : = inf { r, 


r] - x 2 + ^ in (x)(x 2 - 1 ) + t' 2 in (x)(x 2 - 1 ) = cr 0 (x), 
?7 Gl, do G E[x], uf 1 , vlff G E[x], 
deg(u 0 ) < 6 ,deg((^ in ) < 6,deg{{v% in ) < 6 


It is easy to see that taking for all x G M 2 , iq m (x) = 1 and for all x G M 2 , (x) = 0 leads 

to rj — x 2 + zz” in (x)(x 2 — 1) + t' 2 in (x)(x 2 — 1) = 77 — 1 = cro(x) . T hus for q = 1 and for all 
x G M 2 , <7o(x) = 0, we obtain X™ n (q\) < 1. We will see at Prop, f.l, 


since X in \q{) = supjx 2 | (xi,X 2 ) G [—1,1] x [—1,1]} = 1, we conclude that 1 < A in/V (( 7 i) and 

x in7Z ( gi ) = 1. 


that A mt < T X 




Thus, 


IK, 


Finally, we define the relaxed functional F ^ for all w G F 


for all p G P as follows: 


(F ll (w)){p) = sup <j sup {F 1 i i (w)){p),X m ' IZ {p) 
[ieX 


(13) 


As we followed the construction proposed in Section 4.1, the relaxed functional F ^ provides a safe 
over-approximation of the abstract semantics F z . 
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Proposition 4.1 (Safety). The following statements hold: 

1. A int < F X in1Z ; 

2. For all i G T, for all w G F Ff(w) <jr Fj i {w); 

3. For all w G T F$(w) <jr F ll {w). 

An important property that we will use to prove some results on policy iteration algorithm is 
the nronotonicity of the relaxed functional. 

Proposition 4.2 (Monotonicity). 1. For all i G I, w H > Fj^(w) is monotone on F 

2. The function w H > F^(w) is monotone on F (P,Mj . 

Proof. Let p G P. For d £ J^P,lJ, A £ F (P,M+), p G X[x] ni , 7 G S[x] n ° and 7 G M such 
that deg ((p,r l ) + ( 7 , r 0 )) < 2m deg T*, we define the polynomial in x, ,ri {v) := p — p oT l — 

E q epA {q)(v(q) — q) + (p,F) + ('IF 0 )- We define for v G F (p,r) the set R(v) = {(A, p, 7 , 77 ) | 

ip X ’ fl ’' y ’ ri (v) G X[x], A G F (P, M+) , p G T,[x] ni , 7 G X[x] n °, 77 G M}. Now, let us take w,w' € F (p, 
such that w <jr w'. We have: 

= p-poT l - E qS p A (q) (■w(q ) - q) + (p, r l ) + ( 7 , r°) 

= p-p or- E qS P A(q) (w(g) - u/(g) + u/(g) - q) + (p, r*) + ( 7 , r°) 

= p — poT 1 — E qeP \(q)(w'(q) - ?) + (/i,P) + ( 7 ,r°) - E qe pA(g)(«;(?) - w'(q)) 

= + E qe p A (q)(w'(q) - 10(g)) . 

Then, from w <jr w' and the fact that \(q) are nonnegative scalars, if is an SOS 

polynomial, so is as a sum of a SOS polynomial and a nonnegative scalar. Hence, 

we have R(w') C R(w). Finally, we recall that if A C B, then inf# < inf^. We conclude that 
(F- 1 H)(p) < {FF(w'))(p). 

2. The mapping F ^ is monotone as supremum of monotone maps. □ 


From the third assertion of Prop. 


Prop. 3.1 


4.1 if w satisfies F^(w) <jr w then F^(w) <jr w and from 
w* is an inductive invariant and thus 91 C w*. This result is formulated as the following 


corollary. 

Corollary 4.1 (Over-approximation). For all w G F (V, such that F^{w) <j- w then 91 C w*. 


5 Policy Iteration in Polynomial Templates Abstract Domains 

We are interested in computing the least fixpoint 91^ of F 7 ' , 91^ being an over-approximation of 
91 (least fixpoint of F). As for the definition of 91, it can be reformulated using Tarski’s theorem 
as the minimal post-fixpoint: 

= min{wj G F (V,]R) \F n (w) <jr w} . 

The idea behind policy iteration is to over-approximate 9! R ' using successive iterations which are 
composed of 
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• the computation of polynomial template bounds using linear programming, 

• the determination of new policies using SOS programming, 

until a fixpoint is reached. Policy iteration navigates in the set of post-fixpoints of F ^ and needs 
to start from a post-fixpoint w° known a-priori. It acts like a narrowing operator and can be 
interrupted at any time. For further information on policy iteration, the interested reader can 
consult (GGG+051 IGGTZ07| . 


5.1 Policies 


Policy iteration can be used to compute a fixpoint of a monotone self-map defined as an infimum 
of a family of affine monotone self-maps. In this paper, we propose to design a policy iteration 
algorithm to compute a fixpoint of F In this subsection, we give the formal definition of policies 
in the context of polynomial templates and define the family of affine monotone self-maps. We do 
not apply the concept of policies on F ^ but on the functions FS exploiting the fact that for all 
i G 1, FS is the optimal value of a minimization problem. 

Policy iteration needs a selection property , that is, when an element w G F (P,M^ is given, 
there exists a policy which achieves the infimum. In our context, since we apply the concept of 
policies to Fj it means that the minimization problem involved in the computation of FS has 


an optimal solution. In our case, for w G F (V, and p £ P, an optimal solution is a vector 
(A, < 7 ,//, 7) G F (P,R + ) x S[x] x S[x] ni x S[x] n ° such that, using (12), we obtain: 


{Ff[w))(p) = poT 1 + J2 x (q)(™(q) -q)~ (hf 1 ) - + ^ 

q £p ■ (14) 

and deg(u) < 2mdegT*, deg ((n,r l ) + ( 7 , r 0 )) < 2m deg T* 


Observe that in Eq. (14), (F^(w))(p) is a scalar whereas the right-hand-side is a polynomial. The 
equality in this equation means that this polynomial is a constant polynomial. Then we introduce 
the set of feasible solutions for the SOS problem {Fj l (w)){p)\ 


Sol(u;, i,p) = {(A, < 7 , p, 7) G F (P, M+) x E[.x] x E[x] ni x E[x] n ° | Eq. (14) holds} 


(15) 


Since policy iteration algorithm can be stopped at any step and still provides a sound over¬ 
approximation, we stop the iteration when Sol (w,i,p) = 0. Now, we are interested in the elements 
w G F (P,M) such that Sol (w,i,p) is non-empty: 

FS (p,l) = {w Gj(p,l) I Vi el, Mp e P, Sol(u;,*,p) + 0} • (16) 

The notation FS (P,M^ was introduced in }AGG12| to define the elements w G F (V,M^ satisfying 
Sol (w,i,p) / 0. In |AGG12l Section 4.3], we could ensure that Sol (w,i,p) / 0 using Slater’s con¬ 
straint qualification condition. In the current nonlinear setting, we cannot use the same condition, 
which yields a more complicated definition for FS (P, . 

Finally, we can define a policy as a map which selects, for all w G FS (P, , for all i G F and 

for all p G P a vector of Sol (w,i,p). More formally, we have the following definition: 
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Definition 5.1 (Policies in the policy iteration SOS based setting). A policy is a map 7 r : FS (P, MJ 1 —> 

((X xP)4 X (P,M + ) x E[x] x X[x] ni x Sfrc]" -0 ) such that: Vw £ FS ^P,M^, Vi £ X, Vp £ P, 

7 T (w)(i,p) £ Sol(u;,z,p). 

We denote by n the set of policies. For 7 r £ II, let us define tt\ as the map from FS (P, to 

(X x P) £ J(P,M + ) which associates with w £ FS (P,]R) and (i,p) £ X x P the first element of 
Ti{w){i,p) i.e. if ir(w)(i,p) = (A,a,//, 7 ) then 7 T\(w)(i,p) = A. The equality ii\(w)(i,p) = A means 
that when we perform the policy iterations algorithm, we select the vector of Lagrange multipliers 
A associated with the constraints of the form g(x) < w(q ). The purpose of this selection is to 
update the value of vj using the direction A. The other coordinates composing ir(w)(i,p) that is 
c 7 , p, 7 do not serve the policy iterations algorithm but are only used to take in consideration the 
sets X 1 and X° in the computation of Fj z (w)(p). 

As said before, policy iteration exploits the linearity of maps when a policy is fixed. We have 
to define the affine maps we will use in a policy iteration step. With 7 r £ II, w £ FS (P, , * £ X 

and p £ P and A = 7 TA(ie)(i,p), let us define the map ip : X" (P, M as follows: 

v ^ X (<l)v{q) + {F^(w)){p) - Kq)w{q) ■ (17) 

geP geP 


Then, for 7r £ II, we define for all w £ FS (V, , the map <f>™^ from X 1-7 X (V, R^. Let 

v £ X (P, and p £ P: 

K {w \ v )(p) = sup |sup^^ p (u), X in 7 e (p)| . (18) 


Example 5.1. Let us consider Example f.2 and the function ur(gi) = wr(g 2 ) = 2.1391 and 
w°(p) = 0. Then there exists two SOS polynomials JI and a such that, for all x £ M d ; 


(F^(w)){qi) = (0.687xi + 0.558x 2 - 0.0001xix 2 ) 2 +A(gi)(2.1391 -x?) 

+A(g 2 )(2.1391 - x 2 ) - A(p)p(x) - p(x)(1 - x 2 ) + a(x) 
= 1.5503, 


with A(gi) = A(g 2 ) = 0 and A(p) = 2.0331. It means that A, p and a are computed such that 
(0.687x1 + 0.558x 2 - 0.0001xix 2 ) 2 + A(gi)(2.1391 - x\) + A(g 2 )(2.1391 - x 2 2 ) - A(p)p(x) - /Z(x)(l - 
xf) + u(x) is actually a constant polynomial. 

Then (A, p, a) £ Sol(ix°, 1, q\) and we can define a policy -k(w°) such that 7 r(ix°)(l, gi) = (A, p, a) 
and thus 7TA(tc°)(l, gi) = (0, 0, 2.0331). We can thus define for v £ X (P mn ,M), the affine mapping: 
4i gi (u) = ^(gi)u(gi)+A(g 2 )u(g 2 )+A(p)u(p)+(X] 7 ? '('u;))(gi)—A(gi)u)(gi)—A(g 2 )u>(g 2 )—A(p)u;(p) = 
2.1391u(p) + 1.5503. 

Let us denote by X (P, M) the set of finite valued function on P i.e g £ X (P, M) iff g(p) £ M for 
all p £ P. 

Proposition 5.1 (Properties of Pi )W)P )~ Let n £ II, w £ FS (P,Rj and ( i,p ) £ X x P. Let us write 
A = ir\(w)(i,p). The following properties are true: 
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i- Vw,i, P is a ffi ne on F ( F > K ); 

2. iPw,i,p *' s monotone on F ^P, IK) ; 

3. V«gj(p,r), F^{v){p) < <p*' itP (v); 

4 • ^ w (™) = ^W(p)' 

Proof. Let w E FS (P, , i E X, p E P and 7 r E LL 

1. The fact that is affine follows readily from the definition (Eq. ©)■ 

2. The monotonicity of follows from the nonnegativity of n\(w)(i,p). 

3. Let v E F (V, . Since w E FS (V, , there exists (A, a, p,j) E F (P, M+) x E[x] x E[x] ni x 

E[x] n ° such that deg(cr) < 2?ndegT*, deg ((p,F) + ( 7 ,r 0 )) < 2mdegT*: 

{F^{w))(p) =poT l + J2 Mq)(w(q) -q)- (i*,F) - (tf°) + &■ 

geP 

Writing A = TT\(w)(i,p), we get: 

pi,i, P (v) = J2 x (q) y (q) - J2 x (q) w (q) +p°t 1 + J2 x (q)( w (q) - q) 

geP geP geP 

-(p.,r l ) - ( 7 ,r°) + cr 

= poT l + J2 X {q)(v{q) ~q)~ (h, F) - ( 7 , r°) + cr . 

geP 


Finally, 

( pi 1 i,pi v ) ~ p 0 T% ~ 55 x (q)( y (q) -q) + (m, F) + ( 7 , r°) =a, 

geP 


and recall that (Eq. @) 

(FF( w )){p)= inf n 

A,cr,/z,7,77 


S. t. 


rj-poT 1 -J2 x {q)(w(q) -q) + </L F) + ( 7 , r°) = a 

qeP 

A E F (P, M+) , cr E E[x], p E E[x] n % 7 E E[x] n °, rj E M, 
deg(cr) < 2m deg T*, deg((^t,r*) + ( 7 , 7 -°)) < 2m deg T*. 


(19) 


From Eq. (19), (A, a, p, 7 , Tw,i,p{ v )) is a feasible solution of the latter minimization problem and we 

conclude that (FF(v))(p) < %v X }^ hP \ v ) ■ 

4. . „ 

t Pli^ KhP) ( w ) = J2 X (q'* w (q) + C^MK V ) - 55 = (F^(w)){p) . 

geP gSP 

□ 


The properties presented in Prop. 5.1 imply some useful properties for the maps <Iv • 

Proposition 5.2 (Properties of ). Let it E II and w E FS (P,M^. The following properties 
are true: 
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1. is monotone on F (P,IR) ; 

2. F n <? $l {w) ; 

3. K {w) (w) = F n (w) ; 


4■ Suppose that the least fixpoint of is L e J 7 (P, M). 

unique optimal solution of the linear program: 


Then L can he computed as the 


inf 


I p'eP 


v(p) | V (i,p) el x 


<Pi,W0 XtJ> \ V ) ^ v (p)’ y <i G 




, x m (q) < v(q) 


( 20 ) 


LP problem (20) corresponds exactly to the linear program presented in the case of quadratic 
templates |AGG121 Eq. 4.4], 


Proof. Let n £ II and w e FS (P, . 

1. The map ^ is monotone as the map monotone for al Mel and for all p e P, 

and the the fact that the point-wise supremum of monotone maps is also monotone. 

2. Let v e F ^P, and let p e P. Recall that: 


( F (v))(p) = sup i sup (F- i (v))(p), X 

i£l 


mn (p) 


and from the third assertion of Prop. 


5.1 


we have for all i el, F^(v)(p) < by taking the 

supremum over 1 and then the supremum with X m ^(p), we obtain that F^(v)(p) < ^w W \v)(p), 
yielding the desired result. 

3. This result follows readily from the fourth assertion of Prop. 5.1 and the definition of 
(Eq. @). 

4. By Tarski’s theorem and as is monotone, has a least fixpoint in F (V, Rj. Let L 

be this least fixpoint supposed to be finite valued. Now, from Tarski’s theorem and the definition 
of &w^ w \ we have: 


L = infju | <jr u} 

= inf jv | M (i,p) el x P, 


<plwv )M ( v ) < v (p )> e 


■ A" 


r\lZ, 


(■ Q ) < «(?)} • 


Let us suppose that there exists a feasible solution v such that Y^q£F v (q) < pL(q). Note that 


since X m ^ <jr v, ^2 qe pV is finite. Then we have inf{u,L} < L and inf{u,L} ^ L. As is 


7r(w) 


monotone and as v and L are feasible, we have <f> 


t(w) 


minimality of L. We conclude that L is the optimal solution of Linear Program (20). 


(inf{u,L}) < inf{u,L}. This contradicts the 

□ 


Remark 1. We recall that the linear constraints in Problem (20) come from the use of the function 
defined at Equation © which is affine on the variable v. The linear forms are defined from 
the vector of Lagrange multipliers A found when we solve the minimization problem involved in 
Equation (12). If we had allowed a vector of SOS polynomials A as vector of Lagrange multipliers, 
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we would obtain a set of polynomial inequalities that we would solve using SOS programming. The 
resulted problem would not have a feasible solution. 

For example, let us consider an SOS polynomial template p, an SOS (non scalar) polynomial A 


and a scalar c. Then, in this case, an analog of Problem (12) would be: 


min{ti(p) G M | A (x)v(p) + c < v{p),M x G M, v(p) > X m (p)} 


We assumed that p is a SOS polynomial template, implying that X m (p) is strictly positive. Since 
X(x) is a non scalar SOS polynomial and v(p) > 0, then v(p)( 1 — A(x)) — c is negative for some 
sufficiently large x. This proves the infeasibility of the problem. 

Recall that a function g : >->• M is upper-semicontinuous at x iff for all (x„) n& ^ converging to 

x, then lim sup n _>. +0O g(x n ) < g(x). 

Proposition 5.3. Let p G P. Then w i— >• F^(w)(p) is upper-semicontinuous on XS ^P,M^ n 
X (P, R). 


Proof. Let 7r G II, w G XS (V, IR) 
elements of X (P, 


fl X (P, M) and p G P. Let i G I. Let (w n )ne n be a sequence of 
to w. Let A = 7 T\(w)(i,p). Since ^ is affine on X (P,M), then 

and finally v >->■ &w^ w \v)(p) is continuous on X (P,R) as a finite 
supremum of continuous functions on X (P,M). Then from the second point of Prop. 5.2 for all 
n G N, F^(w n )(p) < ^w^ w \w n )(p). By taking the limsup, we obtain: lim sup n _>. +00 F ri (w n )(p) < 
limsup n _ >+00 (w n ){p) = $l {w \w){p) = F n (v)(p). □ 


<Pi wp is continuous on X 


5.2 Policy Iteration 

Next, we describe the policy iteration algorithm. We suppose that we have a post-fixpoint w° of 
F n in ^(P,R). 

We detail step by step the algorithm presented in Figure [4] At Line[lJ the algorithm is initialized 
and thus k = 0. At Line|iJ we compute F 1z (w k ) using Eq. ( [13] ) and solve the SOS problem involved 
in Eq. (12). At Line pi if for all i G X and for p G P, the SOS problem involved in Eq. (12) has 


an optimal solution, then a policy 7r is available and we can choose any optimal solution of SOS 
problem involved in Eq. (12) as policy. If an optimal solution does not exist then the algorithm 
stops and return w k 

Then 


can 


Now, if a policy n has been defined, the algorithm goes to Line 
define ^ following Eq. (18[). 


12 


and we 


on templates w k+1 as the smallest fixpoint of K Finally, at Line 


we solve LP problem (|20|) and define the new bound 

k\ 


13 


k is incremented. 


w 


If for some k G N, w k XS ^P, and w k 1 G XS ^P, then the algorithm stops and returns 
. Hence, we set for all l > k, w l = w k . 


Theorem 5.1 (Convergence result of the algorithm presented in Figure [4]). The following state¬ 
ments hold: 


1. For all k G N, w k G X (P, M) and F n (w k ) < w k ; 

2. The sequence (w k )k >o generated by Algorithm^is decreasing and converges; 
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input 


e F ( 


), a post-fixpoint of F n 


1 

2 

3 

4 

5 

6 

7 

8 
9 

10 

11 

12 

13 

14 


output: a fixpoint w = F n (w) ifVfc £N, w k £ FS l 
k= 0 ; 

while fixpoint not reached do 

begin compute the next policy n for the current iterate w k 
Compute F n (w k ) using Eq. (131 and Eq. (12 1 ; 


or a post-fixpoint otherwise 


if 


•jo k £ FS (1 
Define n(w k ) ; 


then 


else 


j return u> fc ; 

end 
end 

begin compute the next iterate w k+1 

Define ^ and compute the least fixpoint w k+1 of from Problem (20 1 ; 

k=k+l; 

end 


is end 


Figure 4: SOS-based policy iteration algorithm for PPS programs. 


3. Let w°° = limfc _>. +00 w k , then < w°°. Furthermore, if for all k £ N, w k £ JLS ^P, 

and if w°° € FS (V, then F 7Z (w°°) = w°°. 

Proof. 1. We reason by induction. We have F^(vP) < w° and w° £ F (P, M) by assumption. Now 
suppose that for some k £ N, F n (w k ) < w k and w k £ F (P, M). If w k ^ FS then w l = w k 

for all l > k and then we have proved the result. Now suppose that w k £ FS ^P, and let us 
take 7 T £ II such that \w k ) = F n (w k ). From induction property w \w k ) < w k and thus 

w k is a post-fixpoint of belonging to F (P,M). Since every post-fixpoint of is greater 

than X m ^' then least fixpoint of d )7rl } 7,; ^ is finite valued and thus it is the optimal solution w k+1 of 
Problem (20). Moreover from the second point of Prop. 5.2, F^(w k+1 ) <jr d* 7 ’^ ' > (ic fc+1 ) and since 


w k+1 is the least fixpoint of \u> k+1 ) then F^(w k+1 ) < w k+1 . This completes the proof and 
for all k £ N, w k £ F (P, M) and F^-(w k ) < w k . 

2. Let k £ N. If w k FS (P, ilj then w k+1 = w k < w k . Now suppose that w k £ FS (P, ]R) and 

let 7T £ II such that \w k ) = F^fw k ), then from the third point of Prop. 5.2, \w k ) = 

F^(w k ) < w k \ the inequality results from the first assertion. Then w k is a post-fixpoint of d *^? 77 \ 

Since w k+1 is the least fixpoint of and ^ is monotone then from Tarski’s theorem 

w k+1 < w k . From the first point, for all k £ N, w k £ J 7 (P, M). Moreover by definition of F 

X in ^ < F n (w k ) for all k £ N, then from the first point (w k ) k>o is lower bounded then it converges 
to some w°°. 

3. If for some k, w k ^ FS and w k ~ 1 £ FS ^P, then w°° = w k and we have 

F R '(w oc ) < w°° from the first point. Now suppose that for all k £ N, w k £ FS (V, m). Since 
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F ^ is monotone then for all k E N, F^'(w°°) < F^(w k ) < from the first point. Now taking 
the limit of the right-hand side, we get F(w°°) < w °°. Now, let k E N and let tt E II such 
that <!> \w k ) = F^(w k ). From the second point, w k+1 < w k and from the monotonicity of 
\ we have w k+1 = \w k+1 ) < \w k ) = F 1z (w k ). By taking the limsup on k , we 

get w°° < limsup fc _^ +00 F n {w k ). As F n is upper-semicontinuous on FS ^P, n J 7 (P, R) then, if 
w°° E FS (V, R^, w°° < limsup fc _ >+00 F n (w k ) < F^(w°°) and so w°° = F' R '(w°°). □ 


5.3 Initialization and templates choice 

In Section [3j we have made the assumption that the template basis was given by an oracle. More¬ 
over, in Algorithm [ 4 J we suppose that we have a post-fixpoint w° E F (P,M) of F Now, we give 
details about the templates basis choice and the computation of a post-fixpoint w° > w°°. The 
templates basis choice relies on the computation of a template basis composed of one element. This 
single template is constructed by the method developed in jAGMlfil and is then completed using 
the strategy proposed in |AGM15I Ex. 9]. The single template computation also permits us to 
compute w°. Actually, the method developed in [ AGM15 ] is constructed by using the definition 
of being a post-fixpoint of F 7? . Indeed, suppose that the templates basis is constituted of one 
template p then w° is a post-fixpoint F 72 if and only if F' l ^(w°)(p) < w°(p). This is equivalent to: 


X inK = inf {r/ | v - P + E v t r t G v ” G s N nin l < , 

3 = 1 


and for all i E X: 


(F^(u; 0 ))(p) = inf rj 


< 


S. t. 


r) — p o T l — A l (w° — p) + (p, r l ) + ( 7 , r°) E S[x] 
A* > 0, p E Fj[x] ni , 7 E E[x] n °, p E M 


By definition of the infimum, it is equivalent to the existence of E S[x] and for all i E X of 
A* > 0, p 1 E E[x] ni , 7 * E E[x] n ° such that: 


w° -p+ vfrf € S[s] 
w° — p o T l — A l (w° — p) + (p, r l ) + ( 7 , r°) E £[x] 


( 21 ) 


Now to find a template, it suffices to find p such that Eq. (21) holds. However, the following two 
issues remain. 


First, without an objective function, p = 0 is a solution of Eq. (21). A workaround to avoid 
this trivial solution consists of optimizing a certain objective function under the constraints given 
in Eq. (21). In }AGM15| . a similar optimization procedure (Problem (13) of |AGM15] ) is used to 
prove a property of the form 91 C {x E W l \ n(x) < a}, for a given real-valued function k. Here, we 
are interested in proving the boundedness of the reachable value set, which corresponds to minimize 
a with k = || • 11 2 • 

Second, finding A* and p satisfying Eq. (21) boils down to solving a bilinear SOS problem, which 
is not easy to handle in practice. Thus, we fix A* = 1 as in Lyapunov equations. We also take 
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w° = 0 since p has a constant part. Finally, to obtain a template p, we solve the following SOS 
problem: 


inf 

peK[3"]2m,«'SlR 

s.t. 


w 


- P = &0 - G ' r 


in 

3'3 


3= 1 


no 


VzGX, p- poT 1 = a 1 -J2 v) r ) - Y l) r \ 


w + p — 


= *l> , 


3 = 1 


3 = 1 


Vj = 1, • • • ,n in , <jj G S[x] , deg(cjjr“) < 2m , 

00 e S[x] , deg(cr 0 ) < 2m , 

Vi G X , ff 1 S E[x] , deg(cP) < 2m.deg T l , 

Vi G X , Vj = 1,... ,rii , fi) G S[x] , deg(/x*-r*-) < 2mdegT* , 

Vi G X , Vj = 1,... ,n 0 , 7 * G £[x] , deg( 7 *r°) < 2mdegX* , 


i/j G S[x] , deg(V') < 2m . 


( 22 ) 


Let (p,w) be a solution of Problem ( |22j ). In [AGM151 Prop. 1], we proved that the set {x G 
| p(x) < 0} defines an inductive invariant. To complete the template basis, we use the strategy 
proposed in |AGM15l Ex. 9], that is, we work with the templates basis {x i-t- xf ,i G [d]} U {p}. 
We thus use the inductive invariant set {p(x) < 0, x? < tc} as initialization i.e. the initial bound 
is w°(q ) = w if q ^ p and w^°\q) = 0 if q = p. As opposed to the approach of |AGM15| . we avoid 
increasing the degree of polynomial p to obtain better bounds on the reachable values set. 


Computational considerations The number of (a-priori unknown) coefficients of the polyno¬ 


mial p (of degree 2m and d variables) appearing in Problem 
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is 


(2 m+d 
d 


Thus, Problem 

/2mdegT ! +d' 


Similarly, the number 
can be reformulated 
SDP variables. Therefore, 
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of coefficients of each cV (resp. /j* and 7]) is ( 2mde s T +d ) 

as an SDP program involving ( 2?? ) / +d ) + Y^ierW + + n o](“ d 

our framework is expected to be tractable when either d or m is small. As mentioned in | ! AG Mini 
Section 4], one could address bigger instances while exploiting sparsity properties of the initial 
system, as in jWKKMOGj . 


6 Experiments 

6.1 Details of the running Example. 

Recall that our running example is given by the following PPS: (X m , X°, (A 1 , A" 2 }, {T^T 2 }), 
where: 

A in = [-1,1] x [-1,1] , f X 1 = {x G M 2 I -x? + 1 < 0} 

X«=R 2 “ d { Y 2 = {x€R 2 |x5-1<0} 
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and the functions relative to the partition {X 1 , A 2 } are: 

T 1 (xi,x 2 ) = 


T 2 (x i,x 2 ) = 


0.687xi + 0.558x2 — O.OOOIX 1 X 2 \ 
-0.292x1 + 0.773x2 J 
and 

0.369x1 + 0.532x2 - O.OOOlx? \ 
-1.27xi + 0.12x2 — 0.0001x1X2^ 


The first step consists in constructing the template basis and compute the template p and bound 
w on the reachable values as a solution of Problem ( |22[ ). We fix the degree of p to 6. The template 
p generated from Matlab is of degree 6 and is equal to 

-1.931348006 + 3.5771x? + 2.0669x^ + 0.7702xix 2 - (2.6284e-4)x? 

—(5.5572e-4)xfx 2 + (3.1872e-4)xix| + O.OOlOxf - 2.4650xf - 0.5073xfx 2 
-2.8032x?x| - 0.5894xixl - 1.4968x1 + (2.7178e-4)xf + (1.2726e-4)xfx 2 
-(3.8372e-4)x?xl + (6.5349e-5)xfx| + (5.7948e-6)xix!f - (6.2558e-4)xl 
+0.5987x? - 0.0168xfx 2 + 1.1066xfx| + 0.3172x?x| + 0.8380x?xl + 0.0635xixl 
+0.4719x1. 

The upper bound w is equal to 2.1343. As suggested in Section +3 we can take the template 
basis P mn = {p, x i->- x 2 ,x e-)- x^}. We write q\ for x ^ x 2 and q 2 for x i->- x|. The basic semi- 
algebraic {x £ l 2 | p(x) < 0, qi(x) < 2.1343, q 2 (x) < 2.1343} is an inductive invariant and the 
corresponding bounds function is w° = (w°(qi), w°(q 2 ), w°(p)) = (2.1343,2.1343,0). 

As in Line [ 4 ] of Algorithm [ 4 J we compute the image of w° by F ^ using SOS (Eq. ©)• We 
found that 


F 1K (w J )(qi) = 1.5503, F K (w u )(q 2 ) = 1.9501 and F K (w")(p) = 0 . 

Since w° G FS Algorithm |ij goes to Line Jbj and the computation of F^(w°) permits to 

determine a new policy The important data is the vector A. For example, for i = 1 and the 

template q\, the vector A is (0, 0, 2.0331). It means that we associate for each template q a weight 
A (q). In the case of A = (0,0,2.0331), X(qi) = 0, X(q 2 ) = 0 and A (p) = 2.0332. For i = 1, the 
template q\ and the bound vector w° , the function ip^ 0 1 fu) = 2.0331u(p) + 1.5503. 

To get the new invariant, Algorithm [l] goes to Line 12 and we compute a bound vector w 1 
solution of Linear Program (20). In this case, it corresponds to the following LP problem: 


We obtain: 


min v(qi) + v(q 2 ) + v(p) 

2.0331ti(p)+1.5503<+qi), 1.0429+p)+1.2235<+q 2 ), 0.9535i;(p)-0.0248<+p), 
0.4578+p)+0.8843<Rgi), 0.2048Rp)+1.9501<u(q2), 0.9985+p)-3.4691e-7<+p), 
l<u(gi), l<+< 72 ), 0<v(p ), (init) 

w 1 (qi) = 1-5503, w 1 (q 2 ) = 1.9501 and w 1 (p) = 0 . 


(i=2) 

0=1) 


Then, we come back to Line [ 4 ] of Algorithm E] and we compute F^iw 1 ) using the SOS program 
Eq. (12). The implemented stopping rule is || F^(w k ) — , u; fc || 0O < le-6 and since H-F^w; 1 ) — 'U! 1 || 00 < 


le-6, Algorithm [4] terminates. The computed sets are presented in Figure i page § Figure [5] 
presents the semialgebraic sets obtained with higher dimensional templates, up to degree 10. Results 
are similar but could lead to different numbers of iterations depending on the degree. In case of 
multiple iterations, the final value is also reached at iterate 1 and is slightly modified by following 
iterations. 
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U_1_!_1_1_!_U 

-1.5 -1 -0.5 0 .5 1 1.5 


The semialgebraic sets denotes templates of degree 7 (blue), 8 (red), 9 (green) and 10 (black). The 
first box denotes the initial bounds as obtained in figure [3] The second one is the one obtained 
after 1 iterations. Except degree 7 that converged in 4 iterations, all others converged in 1. Degree 
10 faced numerical issues and did not allow to refine the bounds without errors. 

Figure 5: Different templates and associated bounds computed with Policy Iterations 


6.2 Benchmarks. 


The presented analysis has been applied to available examples of the control community literature: 
piecewise linear systems, polynomial systems, etc. We gathered the examples matching our criteria: 
discrete systems, possibly piecewise, at most polynomial. In all the considered cases, no common 
quadratic Lyapunov existed. In other words, not only the existing linear abstractions such as 
intervals or polyhedra would fail in computing a non trivial post-fixpoint, but also the existing 
analyses dedicated to digital filters such as |Fer041 IGS071 IAGG12| IR.IGF12] . 

The analysis has been implemented in Matlab and relies on the Mosek SDP solver jAAOOj . 
through the Yalmip jLOlj SOS front-end. Without outstanding performances, all experiments 
are performed within a few seconds per iteration, which makes us believe that a more serious 
implementation would perform better. We recall that the analysis could be interrupted at any 
point, still providing a safe upper bound. 

We next present the examples handled by our SOS policy iteration algorithm: 


Example 6.1. The following example corresponds to ! Fen02\. Ex. 2.1] and represents a piecewise 
linear system with 2 cases handling 3 variables. The initial set is: 


x in = [-i,iy 


The set where the state-variable lies is: 

X° = M 3 . 
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The sets defining the partition of the state-space are: 


X 1 = {(x, y , z) E R 3 |x < 0}, X 2 = {(x, y, z ) E M 3 |x > 0} . 


Finally the dynamics associated to the partition are: 



( x + 0.5 y N 


1 x + Ay + 0.01z\ 

T\x,y,z) = 

—0.3x + 0.8y 
y OAz j 

, T 2 (x, y, z) = 

— 0 . 1 x + 0 . 8 y 
v 0.5 z J 


Example 6.2. We consider the example \Fen02\ Ex. 3.3] which describes a piecewise linear system 
with 4 cases handling 2 variables. The initial set is: 


A in = [—1, l ] 2 . 


The set where the state-variable lies is: 


X° = M 2 . 


The sets defining the partition of the state-space are: 


X 1 = {(x, y) E M 2 |x < -1}, X 2 = {(x,y) E M 2 |x E] - 1,1] Ay > 0}, 

A' 3 = {(x, y) E M 2 |x e] — 1,1] A y < 0}, X 4 = {(x, y) E M 2 |x > 1} . 

Finally the dynamics associated to the partition are: 


T\x,y) = 




0.9x — O.Oly 
O.lx + y + 0.02 


Example 6.3. The following example is the piecewise quadratic system with 2 cases handling 2 
variables lAJl'A Ex. 3]. The initial set is: 

A in = [—1, l ] 2 . 

The set where the state-variable lies is: 

X° = M 2 . 

The sets defining the partition of the state-space are: 


X 1 = {(x, y) E M 2 | — x 4 + x 2 — 1 < 0}, A 2 = {(x, y) E M 2 |x 4 — x 2 + 

Finally the dynamics associated to the partition are: 

ml , \ ( 0.687x + 0.558y — O.OOOlx^A , /o.369x + 0.532y - 

T(I - 9) =I^ -0.292X + 0.773J, )’ T {x ’ y) = ^-1.27x + 0.12j, - 


1 < 0 } 


0 . 0001 x 2 \ 
O.OOOlxy J 
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Example 6.4. The following example is the hand-crafted piecewise polynomial of degree 3 with 2 
cases developed in \AGM15f . The initial set is: 


X in = [0.9, LI] x [0,0.2] 


The set where the state-variable lies is: 


X u = 


The sets defining the partition of the state-space are: 

X 1 = {(x, y ) G M 2 |x 2 + y 2 < 1}, X 2 = {(x, y) G M 2 |x 2 + y 2 > 1} 


Finally the dynamics associated to the partition are: 


r '{x,s)=(^ + ^), T\x,y) = 


' 0.5x 3 + 0.4y 2 N 
—0.6x 2 + 0.3y 2 


The table [l] summarizes the examples considered, the bounds obtained, the degree of the poly¬ 
nomial templates and the number of iterations performed before reaching the fixpoint. 


Examples 


Bounds (ie. x 2 


Degree ff it. 


Running example 2.1 


No good invariant 4 

[1.5503.1.9501] 6 

[1.5503.1.9502] 8 

[1.5500,1.9436] 10 

[1.5503,1.9383] 12 


1 

7 

1 

2 


Example 


6.1 


[3.8260,2.1632,1.0000] 4 1 

[3.7482,1.8503,1.0000] 6 1 

No good invariant 8,10,12 


Example 


6.2 


[1.8359,1.3341] 
[1.5854,1.2574] 
[1.5106,1.2569] 
[1.4813,1.2544] 


4 

6 

8 

10 


2 

5 
4 

6 


Example 


6.3 


[1.5624,1.2396] 4 3 

[1.5581,1.1764] 6 1 

[1.5531,1.1511] 8 1 

No good invariant 10,12 


Example 


6.4 


No good invariant 4,6,8,10 

[1.2100,0.9989] 12 max (10) 


“No good invariant” occurs when the template synthesis fails, i.e. does not provide a sound post- 
fixpoint or some numerical issues occurs during the policy iterations phase. It seems to be due 
to the large size of the SOS problems together with numerical issues related to the interior point 
methods implemented in the relying solvers. 


Table 1: Experiments 
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7 Conclusion 


We proposed an extension of policy iteration algorithms, using Sum-of-Squares programming. This 
extension allows to consider the wider class of disjunctive polynomial programs. In this new setting, 
we showed that we keep the advantage of policy iteration algorithms, while producing a sequence 
of increasingly safe over-approximations of the reachability set. 

As future work, we plan to generalize this algorithm to programs involving non-polynomial 
updates, including square roots, divisions as well as transcendental functions. The computational 
method developed in the paper could be also generalized to other classes of nonlinear switched 
systems involving either random or temporal switching. 
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